Tarea_3

# Librerias a utilizar
library(phyloseq)
library(microbiome)

## Loading required package: ggplot2

## 
## microbiome R package (microbiome.github.com)
##     
## 
## 
##  Copyright (C) 2011-2022 Leo Lahti, 
##     Sudarshan Shetty et al. <microbiome.github.io>

## 
## Attaching package: 'microbiome'

## The following object is masked from 'package:ggplot2':
## 
##     alpha

## The following object is masked from 'package:base':
## 
##     transform

library(ggplot2)
library(vegan)

## Loading required package: permute

## Loading required package: lattice

## 
## Attaching package: 'vegan'

## The following object is masked from 'package:microbiome':
## 
##     diversity

library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

Explorando el objeto phyloseq

data("dietswap", package = "microbiome")
filo_sec <- dietswap
filo_sec

## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 130 taxa and 222 samples ]
## sample_data() Sample Data:       [ 222 samples by 8 sample variables ]
## tax_table()   Taxonomy Table:    [ 130 taxa by 3 taxonomic ranks ]

View(otu_table(filo_sec))

View(sample_data(filo_sec))

View(tax_table(filo_sec))

De acuerdo con el codigo anterior ¿Cuántas muestras y taxones contiene el objeto? - Tiene 222 muestras y 130 taxones

¿Qué variables están disponibles en los metadatos de las muestras? - Son 8, las cuales son: numero de muestra, nombre de la muestra, sexo, nacionalidad, grupo, punto de tiempo, punto de tiempo dentro de grupo, indice de peso.

Curvas de rarefacción

Para realizar las curvas de rarefacción usaremos la función rarecurve de vegan.

Para ejecutar la función es necesario convertir nuestra otu table a un data frame o matriz, y posteriormente invertir la posición de las columnas con la de los renglones.

# generamos un objeto a partir de la otu table
otu_table(filo_sec) -> abundancias

# convertimos el objeto a data frame
data_abundancias <- as.data.frame(abundancias)

# transponemos el data frame para que se invierta la posicion de los renglones con las columnas

t(data_abundancias) -> t_data_abundancias

# Para que la curva de cada sample se distinga la una de la otra le pedi a chat gpt un vector de colores

colores_chat <- c("black", "darkblue", "darkcyan", "darkgoldenrod",
                  "darkgray", "darkgreen", "darkgrey", "darkkhaki",
                  "darkmagenta", "darkolivegreen", "darkorange", "darkorchid",
                  "darkred", "darksalmon", "darkseagreen", "darkslateblue",
                  "darkslategray", "darkslategrey", "darkturquoise", "darkviolet",
                  "midnightblue", "navy","brown","deepskyblue", "skyblue")

# Generamos las graficas utilizando el el data frame transpuesto, variando el rango de renglones (samples) que se estan tomando en cuenta.

rarecurve(t_data_abundancias[1:25, ], 
          step = 25, 
          xlab = "Tamaño muestral", 
          ylab = "Especies",
          tidy = FALSE,
          col = colores_chat,  
          label = FALSE) 

rarecurve(t_data_abundancias[26:50, ], 
          step = 25, 
          xlab = "Tamaño muestral", 
          ylab = "Especies",
          tidy = FALSE,
          col = colores_chat,  
          label = FALSE) 

rarecurve(t_data_abundancias[51:75, ], 
          step = 25, 
          xlab = "Tamaño muestral", 
          ylab = "Especies",
          tidy = FALSE,
          col = colores_chat,  
          label = FALSE) 

rarecurve(t_data_abundancias[76:100, ], 
          step = 25, 
          xlab = "Tamaño muestral", 
          ylab = "Especies",
          tidy = FALSE,
          col = colores_chat,  
          label = FALSE) 

rarecurve(t_data_abundancias[100:125, ], 
          step = 25, 
          xlab = "Tamaño muestral", 
          ylab = "Especies",
          tidy = FALSE,
          col = colores_chat,  
          label = FALSE) 

rarecurve(t_data_abundancias[126:150, ], 
          step = 25, 
          xlab = "Tamaño muestral", 
          ylab = "Especies",
          tidy = FALSE,
          col = colores_chat,  
          label = FALSE) 

rarecurve(t_data_abundancias[151:175, ], 
          step = 25, 
          xlab = "Tamaño muestral", 
          ylab = "Especies",
          tidy = FALSE,
          col = colores_chat,  
          label = FALSE) 

rarecurve(t_data_abundancias[176:200, ], 
          step = 25, 
          xlab = "Tamaño muestral", 
          ylab = "Especies",
          tidy = FALSE,
          col = colores_chat,  
          label = FALSE) 

rarecurve(t_data_abundancias[201:222, ], 
          step = 25, 
          xlab = "Tamaño muestral", 
          ylab = "Especies",
          tidy = FALSE,
          col = colores_chat,  
          label = FALSE)

¿Qué indican estas curvas? - Nos dicen si el numero de muestreos por zona fue suficiente para capturar una riqueza de especies cercana a la realidad.

¿Hay muestras que deberían descartarse por bajo conteo? - Si, pero debido a que realice las graficas para conjuntos de 25 muestras no se sabe exactamente que muestra es, de manera que realice el siguiente codigo para averiguarlo.

# Analizando las graficas podemos notar que en el segmento    que va de 51 a 75 hay una muestra que no llega a su         asintota, de manera que trabajare con ese segmento.

t_data_abundancias[51:75,] -> rango 

# Inicializamos la variable x en 50 para que se sume con i y el contador inicie en 51
x <- 50
for (i in 1:25) {
  sum(rango[i,]) -> tamaño_muestral
  x+i -> muestra
  print(paste0("El tamaño muestral es: ",tamaño_muestral,
                " para la muestra ",muestra))
}

## [1] "El tamaño muestral es: 10819 para la muestra 51"
## [1] "El tamaño muestral es: 10979 para la muestra 52"
## [1] "El tamaño muestral es: 8412 para la muestra 53"
## [1] "El tamaño muestral es: 23529 para la muestra 54"
## [1] "El tamaño muestral es: 16757 para la muestra 55"
## [1] "El tamaño muestral es: 1776 para la muestra 56"
## [1] "El tamaño muestral es: 12207 para la muestra 57"
## [1] "El tamaño muestral es: 10503 para la muestra 58"
## [1] "El tamaño muestral es: 6531 para la muestra 59"
## [1] "El tamaño muestral es: 14325 para la muestra 60"
## [1] "El tamaño muestral es: 13984 para la muestra 61"
## [1] "El tamaño muestral es: 8326 para la muestra 62"
## [1] "El tamaño muestral es: 11732 para la muestra 63"
## [1] "El tamaño muestral es: 9614 para la muestra 64"
## [1] "El tamaño muestral es: 8529 para la muestra 65"
## [1] "El tamaño muestral es: 8869 para la muestra 66"
## [1] "El tamaño muestral es: 10920 para la muestra 67"
## [1] "El tamaño muestral es: 8168 para la muestra 68"
## [1] "El tamaño muestral es: 14956 para la muestra 69"
## [1] "El tamaño muestral es: 9328 para la muestra 70"
## [1] "El tamaño muestral es: 10521 para la muestra 71"
## [1] "El tamaño muestral es: 10391 para la muestra 72"
## [1] "El tamaño muestral es: 10011 para la muestra 73"
## [1] "El tamaño muestral es: 12859 para la muestra 74"
## [1] "El tamaño muestral es: 5841 para la muestra 75"

La muestra que se deberia de descartar es la numero 56, ya que cuenta con 1776 individuos y no llega a su asintota. Esto lo comprobamos en la grafica, ya que es la unica linea cuyo tamaño muestral es menor a 5000.

Diversidad alfa (𝛼)

Calcular y graficar los siguientes indices - Riqueza - Shannon - Simpson

# primero recortareos a aquellas taxa que tienen menos de 1 de abundancia
filosec.prune <- prune_taxa(taxa_sums(filo_sec) > 1, filo_sec)

# Riqueza
plot_richness(filosec.prune, x = "nationality", measures = "Observed", color = "sex")

# La grafica nos muestra la riqueza de acuerdo a la         nacionalidad y ditingue los sexos por colores.

# Shannon
plot_richness(filosec.prune, x = "nationality", measures = "shannon", color = "sex")

plot_richness(filosec.prune, measures="Shannon", x="sex", color="bmi_group")

# Simpson
plot_richness(filosec.prune, x = "nationality", measures = "simpson", color = "sex")

plot_richness(filosec.prune, measures="simpson", x="sex", color="bmi_group")

¿Qué interpretas de estas gráficas? ¿Hay diferencias notorias entre grupos?

Filtrado y transformación

Aplica un filtrado para quedarte solo con los géneros más abundantes de acuerdo a un criterio dado

# Primero crearé un objeto que contenga la riqueza basal y haré comparaciones contra otros objetos
filtro_basal <- prune_taxa(taxa_sums(filo_sec) > 1, filo_sec)
estimate_richness(filtro_basal, measures = "Observed") -> Riqueza_i # Riqueza inicial

filtro_abundancias <- prune_taxa(taxa_sums(filo_sec) > 300, filo_sec)
filtro_abundancias # Este objeto tiene las taxa que tienen una abundancia de al menos 300 en todos los sample

## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 98 taxa and 222 samples ]
## sample_data() Sample Data:       [ 222 samples by 8 sample variables ]
## tax_table()   Taxonomy Table:    [ 98 taxa by 3 taxonomic ranks ]

estimate_richness(filtro_abundancias, measures = "Observed") -> Riqueza_f # Riqueza final

# Cálculo del cambio de la riqueza de acuerdo al filtrado por abundancias
Riqueza_i-Riqueza_f -> cambio # La riqueza de A debe ser más grande que la de B porque esta considerando más taxa

(cambio/Riqueza_i)*100 -> porcentaje_cambio
min(porcentaje_cambio) # Esto quiere decir que en Riqueza_f tego a los taxones que representan al menos

## [1] 4.901961

# un 4.9% - 100% = 95.6% de la riqueza en los sample
max(porcentaje_cambio) 

## [1] 12.61261

# 12.6% - 100% = 87.4%

# De manera que los generos que representan de un 87.4% a 95.6% de la riqueza se encuntran en el objeto filtro_abundacias
# Dentro del cual tenemos taxa que tienen al menos 300 de abundancia en total 
filtro_abundancias

## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 98 taxa and 222 samples ]
## sample_data() Sample Data:       [ 222 samples by 8 sample variables ]
## tax_table()   Taxonomy Table:    [ 98 taxa by 3 taxonomic ranks ]

View(otu_table(filtro_abundancias))

###Diversidad beta Realizar una ordención PCoA utilizando distancia Bray-Curtis

# Primero usamos la funcion distance para obtener las distancias con el metodo Bray-curtis
ditancia_braycurtis <- distance(filo_sec, method = "bray")

# Ahora obtenemos las ordenadas
datos_ordenadas_pcoa <- ordinate(filo_sec, method = "PCoA", distance = ditancia_braycurtis)

# primero hacemos el grafico base generando los puntos y separandolos por color de acuerdo a las variables de los metadatos.
# Aquí mismo se generan los ejes de x y y, seleccionado el punto en el que intersectarán.
# Ademas se añadiran los elipses.

pcoa_plot<- plot_ordination(filo_sec, datos_ordenadas_pcoa, 
                             type = "sample", 
                             color = "group") +
  geom_point(size = 0.5) +  
  geom_hline(yintercept = 0, color = "darkgreen", linetype = 5) +
  geom_vline(xintercept = 0, color = "navy", linetype = 5) +
  stat_ellipse(aes(group = group), level = 0.8, color = "darkred", linetype = "dashed") # Los elipses generados son de acuerdo a los elementos dentro de la variable group, es decir : DI, HE y ED.

# Por ultimo se modificaran las etiquetas de los ejes y el titulo de la grafica
pcoa_plot +                         
  ggtitle("PCoA (Bray-Curtis)") +               
  labs(x = paste0("PCoA1 (", round(datos_ordenadas_pcoa$values$Relative_eig[1] * 100, 1), "%)"),
       y = paste0("PCoA2 (", round(datos_ordenadas_pcoa$values$Relative_eig[2] * 100, 1), "%)")) +  # Etiquetas con porcentaje de varianza
  theme_bw() +                                 
  theme(plot.title = element_text(hjust = 0.5)) 

¿Los grupos se separan visiblemente? - no, hay un area muy grande en la que convergen las 3 elipses generadas.

¿Qué podría estar causando esas diferencias? - La ubicación geometrica que tiene cada muestra dentro de las coordenadas principales del grafico. - En este caso las coordenadas principales no son muy explicativas ya que la PCoA 1 explica tan solo un 42% de la varianza, mientras que la PCoA2 explica el 17.2%, eso implica que los grupos DI, HE y ED no son lo suficientemente explicativos como para realizar una agrupación eficiente de las muestras brindadas.

# calculo de de la abundancia de cada taxon
abundancias_taxones <- taxa_sums(filo_sec)

# Abundancias ordenadas de mayor a menor
abundancias_ordenadas <- sort(abundancias_taxones, decreasing = TRUE)

# gráfica de barras rank-abundance
par(mgp = c(0, 1, -0.5))
barplot(abundancias_ordenadas, 
        main = "Rank-Abundance", 
        xlab = "Taxones", 
        ylab = "Abundancia total", 
        col = "sienna", 
        las = 2,  
        cex.names = 0.5,
        cex.axis = 0.7)

Graficas de abundancia por taxón

# Podemos usar únicamente la funcion plot bar
plot_bar(filo_sec, fill = "Phylum")

# o complementarla
plot_bar(filo_sec, fill = "Phylum") + 
  geom_bar(aes(color=Phylum, fill=Phylum), stat="identity", position="stack")

# separada por grupos
plot_bar(filo_sec, fill = "group") + 
  geom_bar(aes(color=group, fill=group), stat="identity", position="stack")

# Agrupe solo a nivel de phylum porque a nivel de genero o familia se ve un desastre.

¿Hay algún phylum que domine? si, los firmicutes.

¿Se observan diferencias entre grupos? Los HE (azules) tienen una mayor abundancia que los ED (verdes), mientras que los DI (rojo claro) son los que tienen menor abundancia.

Exportar resultados

# Exporta alguna tabla de abundancias o métricas de diversidad a CSV.